#!/bin/bash -e

function info() {
echo Usage: `basename $0` 'in.vcf methy_result fields'
exit 2
}

while getopts  ":p:f:" opt; do
	case  $opt  in
		p) out_prefix=$OPTARG;;
		f) suffix=$OPTARG;;
		*) info;;
	esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi

. $var

# cat wgEncodeHaibMethyl450CpgIslandDetails.txt|grep '^cg'|awk -F',' '{print $12 "\t" $13-1 "\t" $13 "\t" $1 }' |sed '1,8d' > wg450.bed
# cat wg450.bed |sort -k1,1 -k2,2n > wg450.sort.bed
# bgzip -c bed > bed.gz
# tabix -pbed bed.gz

# for i in $(cat cg.id.3722.txt); do cat wg450.sort.bed|grep `echo $i`; done > cg.bed
# cat cg.bed |sort -k1,1 -k2,2n > cg.sort.bed
# bgzip -c bed > bed.gz
# tabix -pbed bed.gz

methyl_id=methyl_id
beta_id=beta
# interval_3722=$data_path/intervals/1/cg.3722.b37.sort.bed.gz
# interval_450=$data_path/intervals/1/wg450_b37.sort.bed.gz
# interval_10k=$data_path/intervals/1/cg.10k.b37.sort.bed.gz
bed450=$data_path/intervals/1/450.clean.$genome_assembly.bed.gz
vcf_header=$data_path/methy/vcf_header
beta=beta.txt
b10k=$data_path/methy/beta/b10k_0.txt
b3k=$data_path/methy/beta/b3k.txt

vcf=$1

cat $2|cut -f$3|sed 1d > $beta


cat $vcf |sed /^##vcfProcessLog/d > $vcf.tmp && mv $vcf $vcf.tmp1 && mv $vcf.tmp $vcf && rm $vcf.tmp1
bgzip -c $vcf > $vcf.gz
tabix -pvcf $vcf.gz

# methylation id annotate
bcftools annotate -R$bed450 -a$bed450 -h$vcf_header -o$out_prefix.cg450.vcf -Ov -c CHROM,FROM,TO,$methyl_id $vcf.gz

# beta annotate
if test -e $beta; then
bgzip -c $beta > $beta.gz
tabix -b2 -e2 $beta.gz

bcftools annotate -a$beta.gz -h$vcf_header -o$out_prefix.beta.vcf -Ov -c CHROM,POS,$beta_id $out_prefix.cg450.vcf

$java_run/snpsift extractField $out_prefix.beta.vcf $methyl_id $beta_id|sed 1d > $out_prefix.cg450.beta.txt


echo -e '\tbeta\ncancer\tna' |cat - $out_prefix.cg450.beta.txt > $out_prefix.cg450.beta.pd.txt

rf_methy_011017.py $b10k $out_prefix.cg450.beta.pd.txt 200 4
rf_methy_011017.py $b3k $out_prefix.cg450.beta.pd.txt 200 4


bk.sh -tproject/methylation/test/`date +%Y-%m` predict_score.txt

fi

. $cmd_done